#if !defined(NR_H)
#define NR_H

#include <memory.h>
#include <fstream>
#include <cstdlib>
#include <iostream>
using namespace std;
#include <cmath>


//void   inverse_jac(double **invin, double **invout, double &det, int ndim);
//void   jacobi (double **a, int n, double *d, double **v);
//void   lubksb (double **a, int n, int *indx, double *b);
//void   ludcmp (double **a, int n, int *indx, double &d);
//void   inverse(double **invin, double **invout, int ndim);
double gasdev (void);
void initiate_ran1(const int l_idum);
double ran1   (void);
void   frprmn(double *p, int n, double ftol, int *iter, double *fret, double (*func)(double *), void (*dfunc)(double *, double *), string& directory);
//void   frprmn(double *p, int n, double ftol, int *iter, double *fret, double (*func)(double *), void (*ofunc)(double *, string &, int), void (*dfunc)(double *, double *), string& directory, double *pstruc, double *pstrucval);
void   dlinmin(double *p, double *xi, int n, double *fret, double (*func)(double *), void (*dfunc)(double *,double *));
void   mnbrak(double *ax,double *bx,double *cx, double *fa,double *fb,double *fc, double (*func)(double));
double dbrent(double ax, double bx, double cx, double (*f)(double), double (*df)(double), double tol, double *xmin);
double df1dim(double x);



double brent(double ax,double bx,double cx,
			 double (*f)(double),double tol, double *xmin);
double f1dim(double);
void linmin(double *p,double *xi,int n,
			double *fret,double (*func)(double*));
void mnbrak(double *ax,double *bx,double *cx,
			double *fa,double *fb,double *fc,
			double (*func)(double x));
void powell(double *p,double **xi,int n,double ftol,
			int *iter,double *fret,
			double (*func)(double*));


template<class T> void newA(T* &p, int nr)
	{
	if(nr==0) nr = 1;
	p = new T[nr];
	p--;
	}

template<class T> void deleteA(T *p)	
{	
  delete[]  (p+1);
}


template<class T> void newM(T** &p, int nr, int nc)
	{
	if(nr==0) nr = 1;
	if(nc==0) nc = 1;

	 p = new T*[nr+1];
	*p = new T [nr*nc];

	p[1] = *p - 1;
	for(int n=2;n<=nr;n++)	p[n] = p[n-1] + nc;
	}

template<class T> void deleteM(T **p)
{	delete[] *p;	delete[]  p;	}


template<class T> void newSymmA(T* &p, int nr)
{
p = nr + new T[2*nr+1];
	
}

template<class T> void deleteSymmA(T* &p, int nr)
{
delete (p-nr);
}


template<class T> bool SameA(T *a1, T *a2, int N)
	{
	if(memcmp(a1+1, a2+1, N*sizeof(T))==0) 
		return (true);
	else
		return (false);
	}


template<class T> bool SameM(T **a1, T **a2, int N1, int N2)
	{
	if(memcmp(*a1, *a2, N1*N2*sizeof(T))==0) 
		return (true);
	else
		return (false);
	}

// copy content of a2 to a1
template<class T> void CopyA(T *a1, T *a2, int N)
{	memcpy(a1+1, a2+1, N*sizeof(T));	}


template<class T> void CopyM(T **a1, T **a2, int N1, int N2)
{	memcpy(*a1, *a2, N1*N2*sizeof(T));	}


class SAFEifstream : public ifstream
{
public : 
	SAFEifstream(void)           { }
	SAFEifstream(char *FileName) { open(FileName); }
   ~SAFEifstream(void)           { }

	void open(char *FileName)
	{
	ifstream::open(FileName);
	if(!ifstream::good())
		{
		cerr << "Error opening file : " << FileName << endl;
		exit(1);		
		}
	}
};


inline double SignOf(double x)
	{
	if(x>0) return ( 1.0);
	else    return (-1.0);
	}

inline double fabs(double x, double &y)
	{
	y = SignOf(x);

	return (fabs(x));
	}

#endif // !defined(NR_H)